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In this paper we consider the problem of bootstrapping a class of 
spatial regression models when the sampling sites are generated by 
a (possibly nonuniform) stochastic design and are irregularly spaced. 
It is shown that the natural extension of the existing block bootstrap 
methods for grid spatial data does not work for irregularly spaced 
spatial data under nonuniform stochastic designs. A variant of the 
blocking mechanism is proposed. It is shown that the proposed block 
bootstrap method provides a valid approximation to the distribu- 
tion of a class of M-estimators of the spatial regression parameters. 
Finite sample properties of the method are investigated through a 
moderately large simulation study and a real data example is given 
to illustrate the methodology. 

1. Introduction. Irregularly spaced spatial data occur frequently in many 
applications and seem to be the rule rather than the exception (cf. [15]). Dis- 
tributional properties of estimators and test statistics based on irregularly 
spaced spatial data are typically influenced by intricate interactions among 
the design density, the spatial asymptotic structure and the spatial cor- 
relation structure of the underlying process, and have complicated forms. 
Resampling methods can play a very useful role in the statistical analysis of 
such spatial data because these methods do not require analytical derivations 
of the exact forms of such interactions and their effects on the (asymptotic) 
distribution of a statistic. Although work on resampling methods for spatial 
data (cf. [13]) began at about the same time as that for temporal processes 
(cf. [7]), much less seems to be known about resampling methods for spatial 
data that are irregularly spaced. In a pioneering work, Politis, Paparoditis 
and Romano [29] developed a subsampling method for irregularly spaced 
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spatial data generated by a homogeneous Poisson process. Politis, Paparo- 
ditis and Romano [30] formulated a version of the spatial block bootstrap, 
also under the same framework. In this paper we consider the problem of 
formulating a unified spatial block bootstrap method for irregularly spaced 
spatial data in a more general framework where the data-sites may have a 
nonuniform concentration across the sampling region and where the data 
generating mechanism may have a nontrivial infill sampling component to 
it. Both these features are important in applications involving irregularly 
spaced spatial data, as (i) the total number of sampling sites and the vol- 
ume of the sampling region may not be of a comparable order of magnitude 
(as required by a homogeneous Poisson process) and (ii) a higher concen- 
tration of sampling sites may be located around the "hot-spots" (leading to 
a nonuniform density of sampling sites) . 

To highlight the peculiarities associated with formulation of a valid block- 
ing mechanism under a nonuniform sampling design, we consider the follow- 
ing example. Let R he a sampling region in the plane. First suppose that 
the sampling sites are generated by a square-grid, for example, the integer 
grid Z^, where Z = {0, ±1, . . .} denotes the set of all integers. For this case, 
Politis and Romano [31] give a formulation of the block bootstrap, where 
(overlapping) blocks are formed by considering translations of a (square- 
)template by points on the grid. Thus, in the regular grid case, one uses the 
data-locations themselves to define the blocks. A random sample of these 
blocks then yields the bootstrap observations. 

Now consider the case where the data-sites are generated by a stochastic 
design and hence are irregularly spaced. In analogy to the regular grid case, 
one may form the blocks by using translations of a given (square-)template 
of a suitable size by the sampling sites themselves [see Figure 1(b)] and 
resample from these blocks to generate the "bootstrap observations." We 
call this version of the spatial block bootstrap the data-site-shifted block 
bootstrap method or the DSSBB method, in short. Although the DSSBB 
method is a "natural" extension of the standard spatial block bootstrap 
method for regularly spaced grid data to the irregularly spaced data case, 
it turns out that the DSSBB method may fail when the spatial sampling 
density is nonuniform. In Section 4 we construct an example of a nonuniform 
spatial sampling density and show that the DSSBB estimator of the variance 
of the sample mean is inconsistent for a wide range of block sizes. As a result, 
the DSSBB is not suitable for irregularly spaced spatial data when the design 
density is nonuniform. 

As a remedy, we also present an alternative formulation of the block- 
ing mechanism and establish its validity for irregularly spaced spatial data, 
allowing nonuniformity of the spatial sampling design and also allowing non- 
stationarity of the underlying spatial process. The alternative formulation 
introduces an extraneous regular grid and defines the blocks by considering 





Fig. 1. (a) Sampling region and the partition R„ = UkeK ^i(k) Wf- (4-3)]. For sim- 
plicity, we write i?„(k) as _Rn(l), . . . , _Rn(12) in the figure. [b(i)] A few representative 
translates of a complete block, say _R„(1), with the lower left points on the sampling sites. 
These translates are used to construct the DSSBB version of the spatial process on the 
complete blocks. [b(ii)] A few representative translates of an incomplete block, say _R„(10), 
with the lower left points on the sampling sites. These translates are used to construct the 
DSSBB version of the spatial process on the incomplete blocks. [c(i)] A few representative 
translates of a complete block, say _R„(1), with the lower left points on the introduced reg- 
ular grid. These translates are used to construct the GBBB version of the spatial process 
on the complete blocks. [c(ii)] A few representative translates of an incomplete block, say 
-Rn(lO), with the lower left points on the introduced regular grid. These translates are used 
to construct the GBBB version of the spatial process on the incomplete blocks. 



translates of a suitable template by points on this grid [see Figure 1(c)]. Al- 
though the alternative blocking mechanism is less natural compared to the 
version based on the data-site-shifted blocks, it circumvents the problems 
associated with the natural blocking mechanism of the DSSBB method. We 
refer to this new version of the block bootstrap method as the grid-based 
block bootstrap method or the GBBB method, in short. 

In this paper we investigate properties of the GBBB method for a class of 
spatial regression models under a stochastic design that is driven by a collec- 
tion of independent and identically distributed (i.i.d.) random vectors with a 
possibly nonuniform density. The spatial asymptotic structure adopted here 
allows "infilling" of any given subregion of the sampling design and thus, 
provides a more flexible framework than the homogeneous Poisson process 
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formulation that is commonly assumed for modeling irregularly spaced spa- 
tial data. The main theoretical results of the paper show that under mild 
regularity conditions, the GBBB method provides a consistent estimator 
of the asymptotic variance of a class of M-estimators of the regression pa- 
rameters. We also show that the GBBB provides a valid approximation to 
the distributions of these M-estimators. Thus, the GBBB provides a unified 
bootstrap method for irregularly spaced spatial data that gives a valid ap- 
proximation, irrespective of nonuniformity of the data-sites, even in presence 
of infill sampling. Finite sample properties of the GBBB method are studied 
in a moderately large simulation study. As an illustration of the proposed 
methodology, we also consider a real data example involving the grazing 
patterns of elks and landscape characteristics in northern Wisconsin, where 
elks are of high conservation interest. 

Block bootstrap methods for time-series data and for grid spatial data 
have been put forward by Hall [13], Carlstein [7], Kiinsch [17], Liu and Singh 
[26], Politis and Romano [31, 32], Zhu and Morgan [41], Zhu and Lahiri [39], 
among others. Subsampling methods for grid spatial data have been de- 
veloped by Possolo [35], Politis and Romano [33], Sherman and Carlstein 
[37], Sherman [36], Garcia-Soidan and Hall [11], Lahiri [20], Lahiri, Kaiser, 
Cressie and Hsu [23] and Zhu, Lahiri and Cressie [40]. Optimal block sizes 
for the spatial resampling method have been given by Nordman and Lahiri 
[27]. Lee and Lahiri [25] employ spatial resampling to develop a least squares 
method for estimation of covariance parameters. Bootstrap and subsampling 
methods for irregularly spaced spatial data have been formulated by Poli- 
tis, Paparoditis and Romano [29, 30]. See Chapter 12 of [21] for more on 
resampling methodology for spatial data. 

The rest of the paper is organized as follows. In Section 2 we introduce 
the stochastic design and the spatial asymptotic framework. In Section 3 
we describe the spatial regression model and establish limit distribution 
theory for a class of M-estimators under the stochastic design of Section 2. 
These results appear to be new and may be of some independent interest. 
Among other things, the results of Section 3 highlight the effects of the 
nonuniform sampling design and of the spatial asymptotic structure on the 
large sample distribution of the M-estimators of the regression parameter. 
In Section 4 we present the DSSBB method and give an example to show its 
inadequacy under nonuniform stochastic designs. In Section 5 we describe 
the GBBB method and establish theoretical results on the validity of the 
GBBB approximation for the distribution of the M-estimators. Results from 
a moderately large simulation study and the real data example are presented 
in Section 6 and Section 7, respectively. For clarity of exposition, proofs of 
the main theoretical results are presented in Sections 8, 9 and 10. 

2. The spatial sampling design. 
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2.1. Sampling regions. Let f/o be an open connected subset of (—1/2, 1/2]'^ 
containing the origin and let Rq be a Borel set satisfying Uq C Rq C Oq, 
where for any set ^ C M'^, A denotes its closure. We regard Rq as a "proto- 
type" of the sampling region R^. Let {A„} be a sequence of positive numbers 
such that n'^/Xn — > as n — > oo for some e > 0. We assume that the sampling 
region i?„ is obtained by "inflating" the set Rq by the scaling factor A.„ (cf. 
[23]), that is, 

(2.1) Rn = XnRo- 

Since the origin is assumed to lie in Rq, the shape of Rn remains the same 
for different values of n. To avoid pathological cases, we assume that for 
any sequence of real numbers {an} with a„ — > 0+ as n ^ oo, the number 
of cubes of the lattice OnZ*^ that intersect both Rq and Rq is 0((a„)~('^~^-*) 
as n ^ oo. The boundary condition on Rq holds for most regions i?„ of 
practical interest, including common convex subsets of R'^, such as spheres, 
ellipsoids, polyhedrons, as well as for many nonconvex star-shaped sets in 
R"^. (Recall that a set A C M*^ is called star-shaped if for any x € A, the line 
segment joining x to the origin lies in A.) The latter class of sets may have 
fairly irregular shape. See, for example, [37] and [20] for more details. 

2.2. The stochastic sampling design. Let f{x) be a probability density 
function on Rq and let {X„} be a sequence of independent and identically 
distributed (i.i.d.) random vectors with probability density function f{x) 
such that {X„} and {Z(s) : s G M'^} are defined on a common probability 
space {n,J-,P) and are independent. We assume that the sampling sites 
si , . . . , s„ are obtained from a realization xi , . . . , x„ of the random vectors 
Xi , . . . , X„,, by the relation 

(2.2) Sj = A„Xj, l<i<n. 

Since X]^ 5 • • • ? take values in Rq , S]^ , . . . , s^^ take values in the entire sam- 
pling region i?„ = XnRo- A similar spatial stochastic design has been used by 
Hall and Patil [15] in the context of nonparametric estimation of the auto- 
covariance function of a spatial process. Note that by the strong law of large 
numbers (SLLN), the number of sampling sites lying in any given subregion 
A of Rn is asymptotically equivalent to nP{XnXi £ A) =n J^-^a •^(^) 
a result, for a density function /(x) that is not a constant function on Rq, 
one may have different expected numbers of sampling sites in two distinct 
subregions Ai and A2 of Rn having the same volume. Thus, this formula- 
tion allows us to model irregularly spaced sampling sites that possibly have 
nonuniform concentration across different parts of the sampling region Rn . 
In contrast, given the sample size n, the standard approach of modeling 
irregularly spaced sampling sites using a homogeneous Poisson point pro- 
cess allows the sampling sites to have only the uniform distribution over 
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the sampUng region. A second important feature of the stochastic samphng 
design is that it allows the sample size n and the volume of the sampling 
region i?„ to grow at different rates. For a positive design density /, when 
the sample size re grows at a rate faster than vol(i?„), the volume of Rn, the 
ratio of the expected number of sampling sites in a given subregion A of Rn 
to vo\{A) tends to infinity. This corresponds to the case of infill sampling 
in the stochastic design case (cf. [8], Chapter 5, [22]) and it is of interest in 
geostatistical and environmental monitoring applications (cf. [8, 23]). Thus, 
both these features of the present stochastic spatial sampling design offer 
some additional flexibility in handling irregularly spaced sampling sites over 
the standard homogeneous Poisson process formulation and are important 
for practical applications. 

In the next section we describe a class of spatial regression models and a 
class of M-estimators of the regression parameter vector under the spatial 
sampling design described above. 

3. M-estimation in spatial regression models. 

3.1. A class of spatial regression models. We consider the spatial regres- 
sion model 

(3.1) y(s) =w(s)'/3 + Z(s), sgM'^, 

where w(s) iIR"^ ^ is a nonrandom weight function, /3 S is the vector 
of regression parameters, {^(s) :s € M'^} is a zero-mean stationary random 
field (r.f.) on W^, p, d G N = {1, 2, . . .} and where we denote the transpose of 
a matrix A by A' . The weight function w(s) in (3.1) may depend on a set 
of covariates xi(s), . . . ,Xq{s) as well as on the spatial location s. Thus, the 
general form of w(s) is 

w(s) =7(2;i(s),...,Xg(s),s) 

for some function 7 : R''"^'^ MP. Next suppose that the y-process is ob- 
served at finitely many locations si, . . . , in the sampling region Rn C M*^, 
generated by the stochastic sampling design of Section 2, and let : M — s- M 
be a Borel-measurable function satisfying 

(3.2) E^{Z{0)) = 0. 

We define an M-estimator of f3 based on {(y(sj), w(sj)')' : z = l,...,n} 
as a measurable solution to the equation (in t G R*') 

n 

(3.3) 5^w(si)^(y(s,)-w(s,)'t)=0. 

i=l 

Note that ipoix) = x, rr G M, gives the least squares estimator of /?. 



BOOTSTRAP FOR SPATIAL REGRESSION 



7 



Distributional properties of M-estimators depend on the population 
characteristics of the underlying spatial process as well as on the spatial 
sampling design in a nontrivial manner. Here we establish large sample dis- 
tributional properties of the M-estimator /3„ under the stochastic design of 
Section 2. 

3.2. Asymptotic distribution of M-estimators. To state the results, we 
need to introduce some notation. For x = (xi, . . . jXd)' S K"', write ||x||i = 
+ • • • + \xd\ and ||x|| = (xf + • • • + x^)^/^, respectively, for the and i"^ 
norms on W^. For any set A C W^, let vol(A) denote the volume (i.e., the 
Lebesgue measure) of A and let |^| denote the size (i.e., the number of ele- 
ments) of A. For y € M, write y+ = max{y, 0}. Let l(-) denote the indicator 
function and let C, C(-) denote generic (nonrandom) positive constants, de- 
pending on their arguments (if any), but not on n. Unless explicitly stated, 
limits in order symbols are taken by letting n — > oo. 

Next we define the strong-mixing coefficient for the r.f. {Z(-)}. Let Tz{T) = 
a{Z{s) : s G T) be the a-field generated by the variables {Z(s) : s G T}, T C W^. 
For any two subsets Ti and T2 of M'', let a(Ti,T2) = sup{|P(j4 r\ B) — 
P{A)P{B)\:A£Tz{Ti), B(^Tz{T2)}, where ci(ri, Ts) = inf{||x - s||i : x G 
Ti, s G T2}- Then the strong-mixing coefficient of the r.f. {Z(-)} is defined 
as 

a(a;6) = sup{Q!(ri, r2) : (i(ri, r2) > o, Ti,T2 G 7^(6)}, 

(3.4) 

a > 0, 6 > 0, 

where 7^(6) is the collection of all finite disjoint unions of cubes in with 
a total volume not exceeding h. Note that the supremum in the definition of 
a(a; h) is taken over sets Ti, T2 that are hounded. This restriction is important 
for d>2, as shown by Bradley [5, 6]. For simplifying the exposition, we 
further assume [cf. [9]) that there exist constants C, ri G (0,oo) and T2 G 
[0, 00) such that 

(3.5) a(a;6) <C-a~^i6""2 for all a>l, 6>1, 

and that T2 = for d= 1. Thus, for d>2, the strength of dependence in- 
creases with the volumes of the sets Ti and T2. We shall specify the exact 
conditions on ti,T2 in the statements below. 

Conditions. 

(C.l) There exists a sequence of nonsingular matrices {A„} such that 

(3.6) A-^[-Bw(A„Xi)w(A„Xi)'](A-^)'^/7 as?i^oo. 
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and for all h£W^, 

(3.7) a;; 

where i7 is a positive definite matrix and where Q{-) is a p x p matrix- valued 
function on W^. 

(C.2) tp is differentiable and its first derivative ip' satisfies a Lipschitz 
condition of order 7 € (0, 1]. 

(C.3) There exist a 5 > and r G N such that: 

(i) ^|V'(^(0))|2^+'5 < 00, ^|'0'(^(O))|2^+'5 < 00. 

(ii) Ti > (2r — l)(2r + 6)d/5, where ri is as in (3.5). 

(iii) 0<r2<^, if d>2. 

(iv) / Q(h)cj^,(x) dx is positive definite, where (T^(x) = Eip{Z{0))ip{Z{x)). 

(v) xo = Eij'{Z{0))^0. 

(C.4) For some oq G [0, 1/8) and for some i 0, 

=sup{||A-iw(s)|| :s e Rn} = 0(minKo,e„(logn)-2A^-'^)/4^i}). 
(C.5) The probability density /(•) of Xi is positive and continuous on 

Rq. 

We briefly comment on the conditions. The first part of condition (C.l) 
[cf. (3.6)] is formulated to add flexibility in the definition of the normaliz- 
ing matrix for the M-estimator /?. One possible choice of A„ is given by 
[£'w(A„xi)w(AnXi)']^/^, in which case we may take H = Ip, the identity 
matrix of order p. However, this choice may not be the most convenient 
for the verification of (3.7), which is a version of the well-known Grenander 
condition in our spatial regression context. See [2, 12, 22] for more details 
about condition (3.7). Condition (C.2) is a set of smoothness conditions on 
the score function ip and is satisfied by many score functions tp, such as the 
function ^po{x) = x, x G R, giving the least squares estimator of (3, and by the 
log-pseudolikelihood function of certain Markov r.f. error processes. For cer- 
tain nondecreasing score functions, the limit distributions of M-estimators 
f3n are derived in [24] using a different method of proof. Conditions (C.3)- 
(C.5) are used for deriving the limit distribution of the M-estimator /3„ 
under the given spatial sampling design and are similar to the condition im- 
posed in [22]. The choice of r in (C.3) will be specified in the statements of 
the results below. For example, for establishing the asymptotic normality of 
the M-estimator (cf. Theorem 1), we used condition (C.3) with r = 1. Note 
that corresponding to a bounded ip, the requirement on the strong-mixing 
coefficient in Theorem 1 reduces to a{a;b) < Ca~^'^^^^b^^^ , a > 1, 6 > 1 for 
d>2, and to a{a; 00) < Ca~^~'^, a> 1 for d = 1, for some arbitrarily small 
e € (0, 1]. The latter is close to the optimal rate for d=l. 



J w(A„x + h)w(A„x)72(x)dx (A;^i)'^Q(h) 



as n 



■ 00, 
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Next define the possible asymptotic covariance matrices of the normalized 

f^n by 

(3.8) Ec = Xo^H-^^,H-^ forcG(0,oo], 

where 

5^00 = / cr^/)(x)(3(x)(ix and T,c = c~^Ha^{0) + T,oo 

(3.9) 

for c € (0, 00). 

Also, let £(W|A') denote the conditional distribution of a random vector 
W given X = (t({X„ : n > 1}), and let g denote the Prohorov metric (cf. 
[4, 28]) on the set of all probability measures on MP. Then we have the 
following result. 

Theorem 1. Suppose that conditions (C.1)-(C.5) hold with r = 1. Also, 
suppose that n/Xf^ —)■ c for some c € (0, 00] . Then 

Xf/^A'^0n-P)^N{O,E,) for almost a// Xi, X2, . . . , (Px), 
that is, 

(3.10) Q{CiXi/^A'^0n-P)\X),N{O,E,))^O as n ^ 00, a.s. (Px), 
where Hc,c€ (0, cxo], are as defined in (3.8) and (3.10). 

Theorem 1 shows that the M-estimator is asymptotically normal for 
almost all realizations of the sampling design random vectors Xi , X2 , . . . , 
with its asymptotic covariance matrix depending on the relative growth rates 
of the sample size n and the volume, |Po|-^n5 the sampling region P„. Let 
c = lim„^oo '^/-^n- When c S (0, 00), the sample size n and the volume of the 
sampling region i?„ grow at the same rate, and this corresponds to the "pure 
increasing domain asymptotic structure" (cf. [8, 22]) in the present setup. On 
the other hand, when c = 00, the sample size grows at a faster rate than the 
volume of Rn and therefore, any given subregion of i?„ of unit volume may 
contain an unbounded number of sampling sites as n — > 00. Thus, for c = 00, 
the sampling region Rn is subjected to infill sampling and we get a "mixed 
increasing domain asymptotic structure" with a nontrivial infill component. 
Theorem 1 establishes asymptotic normality of the M-estimators /3„ under 
both types of spatial asymptotic structures. Furthermore, in view of (3.8)- 
(3.10), Theorem 1 also shows that the infill component leads to a reduction 
in the asymptotic variance of Pn in the mixed the positive definite 

matrix c^^Ha{0) drops out for the "c = 00" case. 
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Remark 3.1. From the proof of Theorem 1, it follows that under the 
regularity conditions of Theorem 1, the estimating equation (3.3) admits a 
sequence of solutions {/3n} such that for some constant C G (0,cxd), 

(3.11) P.\^{\\/3n-m>CXn''^\lognf) = o{l) asn^oo, a.s. (Px), 

where -P.|x(') denotes the conditional probability given X. In particular, by 
the bounded convergence theorem, this sequence of solutions is "consistent" 
in the sense that it converges to (3 in (the unconditional) probability. When 
the solution to (3.3) is unique, the unique solution also satisfies (3.11) and 
therefore is consistent for f3. However, if (3.3) has multiple solutions. Theo- 
rem 1 applies to any sequence that satisfies (3.11). In practice, when 
(3.3) has multiple solutions, one may identify such a sequence by considering 

J /O 

the solution that is closest to a An -consistent estimator (such as the least 
squares estimator) of (3. 

Relation (3.10) shows that the limit distributions of the M-estimators 
depend in a very complicated way on the correlation structure of the error 
process, on the regression function and on the spatial sampling density /. 
In view of the complicated form of the asymptotic variance and its depen- 
dence on the value of the infill parameter c, which is typically unknown in 
practice, it is important to develop resampling methods that automatically 
adjust themselves to the different forms of the sampling distribution of nor- 
malized Pn under various combinations of these factors and produce valid 
approximations in all cases. 

In the next section we describe the natural extension of the standard 
block bootstrap method, the DSSBB, and give an example to show that it 
may fail if the design density is nonuniform. In Section 5 we describe the 
GBBB method and establish its validity for approximating the sampling 
distribution of M-estimators in the full generality of the spatial sampling 
framework of Theorem 1, including nonuniform densities / and both types 
of spatial asymptotic structures. 

4. The DSSBB method and its properties. For simplicity of exposition, 
we describe the DSSBB method only for d-dimensional cubic sampling re- 
gions for a stationary r.f. A description of the method for a sampling region 
of a general shape and for a spatial process satisfying (3.1) can be given by 
routinely modifying the description of the GBBB method given in the next 
section. Hence, for the rest of this section, suppose that Rq = {—1/2, 1/2)'^ 
is the prototype set for the sampling region 

(4.1) R^ = XnRo=[-^,^y, 

and that {y(s):s G M'^} is a stationary r.f. [i.e., /3 = in (3.1)]. 
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4.1. The DSSBB method. Let be a sequence of positive numbers 
such that 

(4.2) + hn/\n ^0 as n ^ oo. 

For simphcity of exposition, in this section we also suppose that r„, = \n/hn is 
an even integer for all n > 1. Here 6„ determines the volumes of the bootstrap 
blocks. Condition (4.2) says that the volumes of the DSSBB blocks grow to 
infinity, but at a rate slower than the volume of the sampling region 
The basic idea behind the formulation of the DSSBB method is to partition 
the sampling region into (i-dimensional cubes of volume 6^ and define 
a version of the spatial process on each such subregion by sampling from a 
suitable collection of data-site-shifted blocks. Specifically, let 

(4.3) Rn= i?n(k) 

kG/C„ 

be the partition of where ii'„(k) = i?„ n [(k + [0, lY)hn], k G Z"*, and 
/Cn = {k e Z'^ : i?n(k) / 0}. For i?„ of (4.1), it is easy to check that /C„ = 
{{ki, kd)' e 1"^ : -r„ < 2A;i < r„ for i = 1, . . . , d}. Next, let b\^^ [i] = Sj + 
[0, l)'^6n, I <i <n, denote the cubes of volume with their "lower left" 
end-points at the sampling sites si, . . . ,s„ and define 

ll^^ = {i:l<i< n, s, + [0, ifbn C Rn}, 

the index set of all data-site-shifted (observed) blocks of volume 6^ that are 
contained in the sampling region i?„ [cf. Figure 1(b) (i)]. Let {jB**(k):k€ 
ICn} be a random sample drawn with replacement from the collection 

{Bli^\i):i ell^^}. Then the DSSBB version of the Y-process over i2„, is 
given by concatenating the observations in the resampled blocks {;B**(k) : k S 
ICn}- More precisely, let yn{A) = {Y{si) : 1 < i < n,Si ^ A} denote the set of 
observations over a set of A C M°'. Then, in this notation, 

(4.4) yn{Rn) = {Y{si):i = l,...,n}, 

the collection of all observations. The DSSBB version of yniRn) is defined 
as 

(4.5) y**iRn)^ u yniK*m- 

kG/C„ 

The DSSBB version of a statistic T„ = tn[yn{Rn)) is given by replacing the 
observations {y{si) : i = 1, . . . ,n} by the resampled values 3^**(i2„). 

When the sampling region is of a general shape, not all subregions -Rn(k) 
in (4.3) are necessarily cubic [cf. Figure 1(a)]. In this case, we may only use 
a part of a resampled cubic block that is congruent to a given subregion 
[cf. Figure 1(b) (ii)] and combine all the resampled observations to define 



12 



S. N. LAHIRI AND J. ZHU 



the DSSBB version of the y-process on ah of (See the description of 
the GBBB method in Section 5 for more details.) Pohtis and Sherman [34] 
estabhsh validity of the DSSBB method for marked point processes, where 
the sampling sites are generated by a collection of weakly dependent random 
vectors, but with the uniform distribution on the sampling region. In the 
next section we give an example in the two-dimensional case to show that 
if the design density is nonuniform, the DSSBB method fails to provide a 
consistent variance estimator even for the simple case where T„ is the sample 
mean. 



4.2. Inconsistency of the DSSBB method under nonuniformity. Suppose 
that {5^(s) : s G M"^} is a stationary r.f. on the plane (i.e., d = 2) and that the 
prototype set Rq is given by Rq = (—1/2, 1/2)^ [cf. (4.1)]. Further, suppose 
that the sample size n and the sequence {A^} satisfy the relation 

(4.6) n = {Xif. 

Note that (4.6) implies a mixed increasing domain spatial asymptotic struc- 
ture with a nontrivial infill component (which corresponds to the case "c = 
oo" in Theorem 1). Next, suppose that the spatial sampling density / is 
given by 

(4.7) fa{xi,X2) = ga{xi)\-i/2,l/2){x2), 

for a G (4, oo) (to be specified later), where ga{xi) is a symmetric function 
that, on the interval (0, oo), is given by 

a/4, 0<xi<a~\ 

linear, a^^ <xi < 2a~^, 

a/[4(a-3)], 2a~^ <xi<l/2, 

0, xi > 1/2. 

For a large, the sampling design puts one-half of its mass on the thin strip 
(— a~^, a~-'^)(— 1/2, 1/2) and has the other one-half mass on the rest of the 
unit square. Next let 

n 
i=l 

denote the sample mean and let y**{Rn) denote the DSSBB version of the 
observed l^-process based on blocks of size bn- Then, following the descrip- 
tion of Section 4.1, the DSSBB version of r„ is given by 

(4.9) T** = average of the resampled observations in y**{Rn)- 

Suppose that we wish to estimate the variance of Tn using the DSSBB 
method. In Proposition 1 below, we show that the scaled variance of Tn, 

(4.10) (72=A^Var.|x(T„), 



(4.8) ga{xi) 
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has a nonrandom limit cj^ (say) for almost all realizations of the spatial 
stochastic design vectors Xi,X2, . . . , where Var.|x denotes the conditional 
variance given X. Define the DSSBB estimator of the scaled variance cj^ of 
Tn or of the limiting parameter o"^ by 

<T2=A$^Var,(rr), 

where T** is the DSSBB version of T„ given by (4.9) and where Var^, de- 
notes the conditional variance given {i^(s) :s G M?} U {Xj :i > 1}. The fol- 
lowing result gives the large sample behavior of o"^. Recall that in this paper, 
{Q,T,P) denotes the underlying probability space and hence, P denotes the 
unconditional probability measure. 

Proposition 1. Suppose that {y(s):s€M^} is a zero-mean bounded 
stationary r.f. with strong-mixing coefficient ay(-;-) [cf. (3.4)] satisfying 

(4.11) aYia;b)<Ciexp{-C2a)b^'^, a>l,6>l, 

for some constants Ci, C2, C3 G (0, cxo). Also suppose that relation (4.6) holds 
and the spatial sampling density is given by (4.7). Then: 

(i) cr^ a.s. (Px), where cr^ = / a{s)ds J /^(s)(is and a{s) = 
Cov(Z(s),Z(0)), seM^. 

(ii) Suppose, in addition, that cr(0) 7^ and (t(s) > for all s € M^. Then 
there exist a S (4, 00) and rj = r]{a) € (0, 1) such that 

(4.12) limsupP((T^ > r?cj^) < 1. 

n— >oo 

A proof of Proposition 1 is given in Section 9. Note that (4.12) shows 
that (T^ is not a consistent estimator of o"^. Thus, the DSSBB method fails 
to provide a valid approximation even for very smooth functionals of the 
sampling distribution (e.g., the variance) of the sample mean under nonuni- 
form spatial stochastic designs. It can be shown that a similar inconsistency 
result continues to hold for the pure increasing domain spatial asymptotic 
structure, 

n/A^ ^ c G (0, 00), 

and also for other cases of mixed increasing domain asymptotic structure 
where n/A^ — > 00 possibly at a different rate. The failure of the DSSBB 
method seems to be an artifact of the interaction between the nonuniform 
design density and of the additional randomness in the data-site-shifted 
blocks induced by the random locations of the sampling sites. In the next 
section we describe the GBBB method, which untangles the interaction of 
these two factors by introducing an extraneous grid to define the blocks. 
We also show that the GBBB method produces valid approximations for 
nonuniform sampling densities under both spatial asymptotic structures 
and under the full generality of the regression model (3.1). 
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5. The GBBB method and its consistency properties. 

5.1. The GBBB method. 

5.1.1. Description of the GBBB method for a stationary random field. 
For simplicity of exposition, in this section we suppose that {y(s) : s G R''} 
is a stationary r.f. [i.e., /? = in model (3.1)]. Let {bn} be a sequence of 
positive numbers satisfying (4.2), that is, 6„ ^ cxo as ?i — > oo but 6„/A„, 
as n ^ CX3. However, unlike in Section 4.1, here we do not require Xn/bn to be 
an integer. Let ICn = {k G Z*^ : k6„ + [0, l)'^bn n Rn ^ 0} denote the minimal 
set of indexes k G Z"' such that the collection {k6n + [0, l)'^bn ■ k G )Cn} gives a 
covering of Rn by disjoint (hyper-)cubes of sides 6„. For k G /C„, let i2n(k) = 
Rn n {kbn + [0, l^bn} denote the part of Rn covered by k&„ + [0, 1)'^6„. Then 
Rn can be expressed as a disjoint union of the sets i?„(k),k G JCn, that is, 

(5.1) Rn= U fi„(k) 

ke/c„ 

[cf. (4.3)]. Note that for a nonrectangular sampling region Rn, the subre- 
gions iin(k)'s lying near the boundary of Rn need not be rectangular, and 
in general, may have different shapes. The proposed spatial block bootstrap 
method defines a bootstrap version of the ^-process on each of the subre- 
gions i?„(k). To that end, let J„ = {i G Z"' : i + [0, 1)"'6„ C Rn} denote the 
index set of all overlapping hypercubes of the form i + [0, l)^6n that are con- 
tained in the sampling region i?„. Also, let {/k : k G /C„,} denote a collection 
of i.i.d. random variables with the discrete uniform distribution on the set 
In, that is, for each k G JCn, 

(5.2) p,(4 = i) = ^, i(zZn, 

where, recall that, for any set A, \A\ denotes its size and P* denotes the 
conditional probability given a{{Z{s) : s G W^} U {Xj : i > 1}) = ^. As a first 
step, we select the cubic blocks Ik + [0, l)'^6n, k G JCn, using the bootstrap 
variables {/k : k G ICn}- Since each subregion i?n(k) in the partition of Rn in 
(5.1) is contained in a cube of volume 6^^, it is possible to inscribe a "copy" of 
Rn(k) in the resampled block Ik + [0, l)'^bn for every i G In- More precisely, 
we define the kth resampled block by 

(5.3) B*ik) = Rn{k) - kbn + /k, k G JCn. 

Note that as i?„(k) C k6„ + [0, 1)'^6„, Rn(k) - kbn + C Ik + [0, l)'^bn. Fur- 
ther, -B*(k) is congruent to i?„(k), since they differ only by a translation. 
Next let Sn = {si, . . . , s„} denote the collection of data-sites and for any set 
A C M*^, let yn{A) = {Y{s) : s G AnSn} denote the collection of observations 
y(sj) over the set A. Then, like the DSSBB method, the GBBB version of 
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the K-process on i2„ is given by pasting together the observations in the 
resampled blocks {B*(k) : k G /C„} [cf. (4.5)]: 

yn{Rn)= u yiKm- 

kG/C„ 

The GBBB version of a statistic T„ based on y{Rn) is given by replacing 
variables in y{Rn) with the resampled observations in y*{Rn). For example, 
for Tn = n~^^^^iY{si), the sample mean, its GBBB version is given by 
T* = the mean of the resampled values in y*{Rn). 

5.1.2. Description of the GBBB method for the spatial regression model. 
Extension of the GBBB method from the case of a stationary r.f. to the 
regression model (3.1) is done as follows. Let {i3*(k) : k € /C„} denote the 
resampled blocks of the GBBB method, as given by (5.3). To define the 
bootstrap version of the ^-process, we first define the bootstrap version of 
the error process Z{-) on Rn and then use the structural relation of the 
regression model (3.1) to construct the GBBB version of the y-process. To 
that end, let 

Z„(sj) = Y{si) - X(si)'/3n, i = 1,. . . ,n, 

denote the residuals, and for A C W^, let Zn{A) = {Z„(s) : s G AnSn} denote 
the collection of the residuals Zn{si) over the set A, where recall that Sn — 
{si, . . . , s„}. With this, we define the bootstrap version of the error process 
Z{-) over the subregion i?n(k) by 

(5.4) z:(i?„(k))^i„(i?;(k)), kG/c„. 

As in the stationary case (cf. Section 5.1.1), the bootstrap version Z*{Rn) 
of the error variables over the entire sampling region is now obtained by 
pasting together the bootstrap versions Z*{Rn{'k-)) for all k e /C„, that is, 
by 

(5.5) Z*{Rn)= U Z:(i?n(k)). 

kG/C„ 

Next note that as the sampling sites si, . . . ,s„ are irregularly spaced, the 
numbers of sampling sites in the blocks i + [0, l)'^bn and j + [0, l)'^6n for any 
two distinct indexes i,j G In are typically different. Let L^^ = |i3^(k) n Sn\ 
denote the number of observations in the resampled block i?*(k). Then the 
bootstrap sample size under the present resampling scheme is given by 

(5.6) N*^Y1 ^k, 

ke/c„ 

which is random and is, in general, different from the original sample size n. 
However, it can be shown that the expected value of N* is asymptotically 
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equivalent to n. Let {s*:i = 1,...,A^*} denote the collection of sampling 
sites in the resampled blocks {i?*(k) : k € /Cn}, that is, 

(5.7) {st:i = l,...,N*}= U [i?:(k)n5„], 

kG/C„ 

and let {Z*(s|) : i = 1, . . . , N*} denote the corresponding listing of the boot- 
strap error variables in the collection Z*{Rn). Then we define the bootstrap 
"observations" y*(-)'s by 

(5.8) y*(s*) = w{s*yPn + Z:{s*), i = l,...,N*. 

Given a random variable T„ = t„(I>„(i?„); /?), where 'Dn{Rn) = {(wn(si)', Y{si))' : 
i = 1, . . . ,n} stands for the data at hand, the (overlapping) GBBB version T* 
of Tn is given by replacing Pn(i?„) with {(w„(s*)', Y*{s*)y -A = 1, . . . , N*} = 
T>'^{Rn) and tn{-) with tN*{-) in the definition of T„: 

(5.9) T* = tN*{Vl{Rn)-Pn)- 

Next we apply the above block bootstrap method to a normalized version 
of the M-estimator: 

(5.10) Ti„ = Ai„ (/?„-/?), 

1 Icy 

where Ai„ = An and A^ is as in condition (C.l). Existing literature 
on bootstrapping M-estimators of regression parameters with independent 
errors and/or time-series errors suggests that there can be more than one 
way of defining the bootstrap version of the normalized M-estimator Ti„ (cf. 
[10, 18]). Here we follow an approach initially put forward by Shorack [38] in 
the context of bootstrapping M-estimators in regression models with i.i.d. 
errors. To define the bootstrap version of Ti„, define the variables ^^^(k; t) = 
EiI*iw(s*)^(^*(s*)-w(s*)'t)l(s*GS;(k)), tGMP,kG/C„. Also, let 



c„(k) = 



N* 

5:w(s*)v(^:(s*))i(s*Gi?:(k)) 

i=l 



k G /C„, where E'* denotes the conditional expectation given Q. Then we de- 
fine the bootstrap version /3* of /3„ as a measurable solution to the equation 
(in t G : 



(5.11) [5;(k;t)-c„(k)] = 0. 

ke/c„ 

The motivation behind this definition is the following. For t G k G 
Kn, setting 5„(k;t) = Er=i[w(si)^(l'(si) - w(si)'t)l(si G iin(k))J, we may 
rewrite the estimating equation (3.3) (defining the M-estimator (3n) as 

(5.12) ^5„(k;t) = 0. 

kG/Cn 
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By (3.2), the expected value of the sum on the left-hand side of (5.12) at 
the true parameter value t = /3 is zero. However, the bootstrap versions 
5'*(k;t) of S'n(k;t) do not automatically satisfy this model restriction (in 
the bootstrap world). Here we introduced the centering values Cn(k) in the 
bootstrap estimating equation (5.11) precisely to ensure that the analog of 
this unbiasedness condition holds at the bootstrap "true" value (3n- 

With the bootstrap version (3* of the M-estimator defined by (5.11), 
we now define the bootstrap version of the normalized M-estimator Ti„ as 

(5.13) T^^ = AM-(3n). 



5.1.3. An example. For an illustration of the main steps leading to the 
GBBB version of T^*„ of the normalized M-estimator Tin, consider model 
(3.1), with the weight function w(s) = 1, s G and ip{x) = x, x € M. Then 
{y (s) : s G R'^} is a stationary r.f. with EY{0) = fi. In this case, j3n = Yn cor- 
responds to the M-estimator of the mean /? = /i, and Tin to the normalized 
sample mean 

Tin = Ai„(y„ - n), 

where Yn = Yl'i=i ^{^i)- The residuals are given by Z„(sj) = ^(sj) — Yn, 
i = 1, . . . ,n. Let 

Bn{i; k) = Rn(M) - k6„ + i, iT„ 

denote the blocks of "type k" [i.e., blocks congruent to i?„(k)] and let 
/in(k) = £;*Er=i yis^)^Si G B*{k))], k e /C„. Thus, for each k G /C„, /i„(k) 
is the average of all "type k" block sums, where each sum extends over the 
sampling sites lying in the ith "type k" block i?n(i; k), i G The centering 
constants c„(k) are now given by c,„(k) = /in(k) — YnE^,L^, k G /C„, where 
recall that denotes the number of observations in the resampled block 
Bn(k). With this, (5.11) can be rewritten as 



N* 

E 

Li=l 



(y*(s:)-t) 



kG/C„ 



C(k) 



0. 



Hence, the bootstrap version of the M-estimator is given by 

P*n = Y:-N*~' ^ c(k), 

kG/C„ 

where Y* = J2iLiy*{s*)/N* denotes the bootstrap sample mean. Further, 
by (5.13), the bootstrap version of Ti„ is given by 

Tin = ^IniPn " Pn) 



(5.14) 



A^n Y* 



N 



kG/C„ 



C(k) - Yn 



Ai„(y„* - fL„ 
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where /i„ = EkeK;„ /in(k)/iV*. 

5.1.4. Variants of the GBBB method. 

Remark 5.1. Let Kin = {k G : k6.„ + [0, l^hn C i?„} denote the in- 
dex set of all cubes k6„ + [0, l)'^fen that are completely contained in i?„ and 
let K,2n = Kn \ ICin denote the index set of all boundary hypercubes. Then 
i?n(k) = k6„ + [0, lYhn for all k € /Ci„, while for k e /C2n, each i?n(k) may 
have a different shape. As a consequence, the -B*(k)'s are of cubic shape 
for all k G /Cin- However, for the boundary regions i?n(k), k € /C2n, the 
collection {i?*(k) :k S /C2n} of resampled blocks may have different shapes, 
depending on the shapes of i?n(k), k G K,2n- 

Remark 5.2. The observation in the previous remark suggests that we 
may define an alternative and simpler version of the GBBB method, by 
restricting attention only to the cubic blocks while resampling. Let Z„ = 
{i G Z*^ : i + [0, l)'^6n C Rn] be as defined earlier and let the i2,„(k)'s be as in 
(5.1). Then, with 0) = i + [0, 1)'^6„, 0) : i G 2:„} is a collection of 

overlapping cubic blocks of sides 6„. For the simpler version of the GBBB, 
we make use of the variable {/k^k G /C„,} of (5.2), but now select the ran- 
dom sample of cubic blocks {Bn{I\^; 0) : k G /C„}, all having the same (cubic) 
shape. The cubic-block GBBB (call it GBBB-CB) version of the error pro- 
cess Z{-) is now defined by concatenating the residuals on the resampled 
blocks: 

Z*{Rn)= U ^(^n(/k;0)), 

kG/Cn 

where Rn = UkGK;„{k+ [0, l)'^&n} is a covering of the sampling region Rn by 
cubes of sides hn as its building blocks. The GBBB-CB version of a random 
variable r„ = t„(I>„(i?„); /?) is given by 

where {s*:i = 1,...,A^*} is the set of data-sites in the resampled blocks 
S„(/k;0), kG/C„, and P;;(i?„) = {(w(s*)', y*(s*))' : i = 1, . . . , iV*} is the 
GBBB-CB version of P(i?„), and where {Y{s*) -.i = 1, . . . ,N*} = y*{Rn) is 
the collection of GBBB-CB observations defined using the resampled error 
variables in Z*{Rn) as in (5.8). Thus the GBBB-CB version of is defined 
as a solution /?* of the equation [cf. (5.11)] 

^ [5:(k;t)-c„(k)]=0, 

kGA^Ti 

where t G M^, S*{k;t) is the sum of all [w(s*)'0(3>*(s*) - w(s*)'t)] over the 
data-sites s* lying in the kth resampled GBBB-CB block Bn{lk',0), and 
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Cn{k) = EieJ„ E"=i w(sj)V'(in(sj))l(sj- G 5n(i;0)), k G ICn- Further, 

the GBBB-CB version of Tin is given by 

Tin = MniPn " Pn), 

where recall that Ai„ = An^^A^. It can be shown, by recasting the arguments 
in the proof of Theorem 2, that the GBBB-CB method provides a valid 
approximation to the distribution of Tin under the same set of regularity 
conditions as Theorem 2 below. 

Remark 5.3. It is evident from the formulation of the GBBB method 
that, like the time-series case (cf. [7]), we may similarly formulate a nonover- 
lapping version of the GBBB method as well. In this case, for each k G ICn, 
we need to select a block at random, with replacement from the collection 
k) : j G Jn} of disjoint blocks of "type k," where Jn = {j G : (j + 
[0, l)'^)bn C Rn} and -B„(j,k) = RnO^) — k-|- Note that a similar modifi- 
cation can also be made to the overlapping GBBB-CB method to formulate 
a nonoverlapping version of the GBBB-CB method. 

5.2. Theoretical validity of the GBBB. Without loss of generality, we 
shall suppose that {X„},j>i, {Z(s):sgM'^} and the bootstrap variables 
{/k : k G ICn}, n>l, are all defined on a common probability space (0, P) . 
Let P*, and Var^. denote the conditional probability, conditional expec- 
tation and conditional variance given Q = a{{X.n : n > 1} U {^(s) : s G M'^}). 
For simplicity of exposition, we suppose that the solution to (3.3) is unique. 
When this assumption does not hold, the result is valid for a suitable se- 
quence of solutions of (3.11), as in Remark 3.1. The main result of this 
section is the following. 

Theorem 2. Suppose that conditions (C.1)-(C.5) hold with r = 3. Also 
suppose that for some 6 G (0, 1), 

(5.15) b-^ + Xi-^bn = o{l) asn^oo 

and that niQn = 0(1) as oo. Then 

sup|P*(ri*„ G^)-P.|x(ri„ G^)HO in P.\^-prohahility a.s. (Px) 

(5.16) 

under both the cases "c G (0,cx))" and "c = oo," where C is the collection of 
all measurable convex sets in W . 

Theorem 2 shows that under some general regularity conditions on the un- 
derlying spatial process, the GBBB method produces valid approximations 
to the distributions of M-estimators under both pure- and mixed-increasing 
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domain asymptotic structures. Furthermore, unlike the DSSBB method, it 
remains vahd for nonuniform spatial sampling designs. The bootstrap ap- 
proximation generated by the GBBB method adapts itself to the various 
forms of the large sample behavior of the normalized M-estimators arising 
from different combinations of various spatial asymptotic structures, spatial 
sampling designs and weight functions, correlation structures of the under- 
lying error process. As a result, the GBBB method serves as a unified, valid 
block bootstrap method for a wide range of irregularly spaced spatial data. 

A second important aspect of Theorem 2 is its validity for a class of 
strongly dependent spatial data. Note that in the presence of the infill com- 
ponent, the variance of the sample sum grows at a faster rate than the sam- 
ple size. As a result, the infill component of the spatial asymptotic structure 
induces conditions of long-range dependence in the data. Although block 
bootstrap methods are known to fail for certain classes of equi-spaced long- 
range dependent time-series data (cf. [19]), the results of this paper show 
that GBBB produces valid approximations even in the presence of long-range 
dependence here. 

Although Theorem 2 remains valid for a wide range of block sizes [cf. 
(5.15)], the accuracy of the GBBB estimator depends on the particular block 
size employed in practice. At this point, we do not know the theoretical opti- 
mal spatial block size. Note that in addition to the geometry of the sampling 
region and the covariance structure of the underlying random process, the 
optimal block length for irregularly spaced spatial data would also depend 
on the sampling density and the spatial asymptotic structure. Quantifying 
the exact form of interaction among these various factors remains a chal- 
lenging problem for future investigation. In the next section we investigate 
finite sample properties of the GBBB method through a simulation study 
and provide some guidelines for choosing reasonable block sizes in practice. 

6. Simulation results. The factors we vary in the simulation consist of 
the sampling region size, the spatial sampling design, the sample size, the 
range of dependence, and the bootstrap block size in the GBBB method. 
We consider = XnRo, where Rq = (—1/2,1/2]^ and the scaling factor 
Xn = 12 or 24. Two types of probability density function / for the spatial 
sampling design are used. One is a uniform distribution over Rq and the 
other a mixture of two bivariate normal distributions, 0.5iV((0, 0)', /2) + 

0. 51V((l/4, l/4)',2/2), truncated outside Rq, where /2 is a 2 x 2 identity 
matrix. Sample sizes are chosen to be n = 100,400,900 and the sampling 
sites si, . . . ,s„ are generated according to /. Now, we consider the class of 
spatial regression models in (3.1). For simplicity, we consider a simple linear 
regression ^{s) = (l,x(s))' and (3 = (/?o,/3i)', where the covariate x(-) = or 

1. The true regression parameters are set at /3o = 25 and /3i = —5. The error 
process Z{-), on the other hand, follows a zero-mean Gaussian process with 
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a spherical variogram that has a unit sill, range r, and no nugget effect. Two 
values of range of dependence are considered: r = 2,4. Thus the simulated 
observations are generated according to Y{s) = /3o + (3ix{s) + Z{s). Again, 
for simplicity, we focus on the identity score function 'tp{x) = x. To carry out 
the GBBB method, we consider bootstrap block sizes 6„ = 2, 4, 6 for A„, = 12 
and bn = 4, 6, 8, 12 for A„ = 24. We also implement an empirical rule for 
choosing the block size [14]. First we apply the GBBB method to obtain the 
variance estimate tp of the /? estimates. Then for / preselected subregions, 
we apply the GBBB method to obtain the variance estimate ipi. The block 
size that minimizes the mean square error I~^^i=i{'4'i — V')^ is denoted by 
6* and the block size selected is 6* scaled by the root of the scaling factor 
An relative to the subregion size. 

For each combination of the factors A.„, /, n, r, we generate S = 500 sam- 
ples. For each sample, we obtain the ordinary least squares estimates /3„ and 
perform M = 1000 resamples for various bn to estimate the variance of Pn- 
Based on the 5 = 500 samples, we estimate the true variance of Thus, we 
can compute the root mean square errors (MSE) of the variance estimates 
(see Tables 1 and 2). For a given sampling region with a scaling factor A„, 
the root MSE tends to be smaller when the sample size n is larger, the range 
of dependence r is smaller, and the bootstrap block size 6„ is about the same 
as the range of dependence. As the scaling factor A„, increases from 12 to 
24, the root MSE becomes even smaller. The simulation results are fairly 
similar for the two types of spatial sampling design. The empirical rule for 
choosing the block size selects the block size that has the smallest root MSE 
in most cases. 

Further, for each sample we construct a 90% confidence interval for (3q 
and Pi- Since we know the true values of Pq and Pi, we can compute the 
nominal coverage probability, based on the coverage of S = 500 confidence 
intervals for each regression parameter (see Tables 3 and 4). For a given 
scaling factor A^, the coverage probability tends to be larger when the sample 
size n is larger, the range of dependence r is smaller, and the bootstrap 
block size bn is about the same as the range of dependence. As the scaling 
factor Xn increases from 12 to 24, the coverage probability becomes higher, 
approaching the asymptotic confidence level of 90%. The simulation results 
are again similar for the two types of spatial sampling design. The empirical 
rule for choosing the block size selects the block size that has the largest 
coverage probability in most cases. 

For the sake of comparison, we implement the DSSBB method for the 
uniform spatial sampling design. Tables 1 and 3 show that the root MSEs 
using the DSSBB method are larger than those using the GBBB method 
when the scaling factor A^ is small and the two methods give comparable 
root MSEs when the scaling factor A„ is large, for the different sample 
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Table 1 

Estimated root MSB of the GBBB variance estimators, with sample sizes 
n = 100, 400, 900, scaling factors An = 12 with block sizes 6„ = 2, 4, 6, and A„ = 24 with 
block Sizes On — 4,6,8,12, and a spherical variogram of ranges 2,4,10, based on 5 = 500 
simulations, each with M — 1000 bootstrap resamples 



n 


r 






^ 71 


12 








A„ =24 






bn 


GBBB 


DSSBB 


bn 


GBBB 


DSSBB 


/3o 




(3o 


/3i 


/3o 


/3i 


/3o 




100 


2 


2* 


0.011 


0.011 


0.018 


0.015 


4* 


0.008 


0.010 


0.009 


0.013 






4 


0.015 


0.015 


0.019 


0.017 


6 


0.009 


0.013 


0.009 


0.014 






6 


0.019 


0.020 


0.022 


0.022 


8 


0.011 


0.016 


0.011 


0.017 
















12 


0.013 


0.018 


0.014 


0.021 




4 


2* 


0.046 


0.013 


0.054 


0.014 


4* 


0.011 


0.011 


0.018 


0.013 






4 


0.042 


0.017 


0.047 


0.016 


6 


0.013 


0.014 


0.017 


0.014 






D 


0.049 


0.022 


0.055 


0.021 


O 
O 


0.015 


0.017 


0.018 


0.016 
















12 


0.019 


0.022 


0.022 


0.020 


400 


2 


2* 


0.008 


0.003 


0.009 


0.003 


4* 


0.003 


0.002 


0.003 


0.003 






4 


0.009 


0.004 


0.009 


0.004 


6 


0.003 


0.003 


0.003 


0.003 






6 


0.012 


0.006 


0.012 


0.006 


8 


0.004 


0.004 


0.004 


0.004 
















12 


0.006 


0.005 


0.005 


0.006 




4 


2* 


0.040 


0.003 


0.043 


0.003 


4* 


0.009 


0.003 


0.010 


0.003 






4 


0.032 


0.004 


0.036 


0.004 


6 


0.009 


0.003 


0.010 


0.003 






6 


0.038 


0.005 


0.042 


0.005 


8 


0.010 


0.004 


0.011 


0.004 
















12 


0.013 


0.005 


0.013 


0.005 


900 


2 


2* 


0.008 


0.001 


0.008 


0.001 


4* 


0.002 


0.001 


0.002 


0.001 






4 


0.008 


0.002 


0.009 


0.002 


6 


0.002 


0.001 


0.002 


0.001 






6 


0.011 


0.002 


0.011 


0.002 


8 


0.002 


0.002 


0.003 


0.002 
















12 


0.003 


0.003 


0.004 


0.002 




4 


2* 


0.042 


0.001 


0.044 


0.001 


4* 


0.008 


0.001 


0.010 


0.001 






4 


0.033 


0.002 


0.036 


0.002 


6 


0.007 


0.002 


0.009 


0.002 






6 


0.039 


0.002 


0.042 


0.002 


8 


0.008 


0.002 


0.010 


0.002 
















12 


0.011 


0.002 


0.013 


0.003 



The design density is a uniform distribution. Also reported are the estimated root MSE of 
tire DSSBB variance estimators. The block size selected most often by the Hall, Horowitz 
and Jing [14] empirical rule is marked by *. 



sizes, ranges of dependence and bootstrap block sizes. However, the coverage 
probability using the DSSBB method is quite a bit lower than that using 
the GBBB method, except for a few cases when A„ = 24 and n = 400,900. 

Finally, we show the boxplots of the S = 500 GBBB variance estimates 
in Figures 2 and 3. In general, there is less variation and skewness when the 
sample size n is larger, the range of dependence r is smaller, and the scaling 
factor An is larger. There is no distinct difference in the boxplots for the two 
types of spatial sampling design. 
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Table 2 



Estimated root MSB of the GBBB 


variance estimators, 


with sample sizes 




n = 100, 400, 900, scaling factors A„ = 12 


with block 


sizes bn = 


2,4,6, and A„ = 24 


with 


block sizes bn 


= 4,6,8,12, anrf 


a spherical variogram 


of ranges 2,4, 10, based on S - 


= 500 




simulations, 


each with M = 1000 bootstrap resamples 








A„ = 12 






A„ = 24 




n r 


bn 


Fir. 

Po 


Pi 


On 


Po 


/3i 


lUU z 


A 


0.014 


0.012 


4 


0.010 


0.011 




4 


0.018 


0.016 


6* 


0.012 


0.014 




6 


0.022 


0.022 


8 


0.014 


U.Ul ( 






— 


— 


12 


0.018 


0.022 


A 


9 


0.057 


0.013 


4 


0.021 


0.011 




4* 


0.053 


0.015 


6* 


0.022 


O.Ui4 




6 


0.062 


0.018 


8 


0.025 


U.UID 










12 


0.031 


0.021 


400 2 


2* 


0.013 


0.003 


4 


0.004 


0.003 




4 


0.014 


0.004 


6* 


0.005 


0.004 




6 


0.017 


0.006 


8 


0.006 


0.005 










12 


0.007 


0.006 


4 


2* 


0.058 


0.003 


4 


0.017 


0.003 




4 


0.050 


0.004 


6* 


0.018 


0.004 




6 


0.057 


0.005 


8 


0.020 


0.004 










12 


0.024 


0.006 


900 2 


2* 


0.012 


0.001 


4* 


0.004 


0.001 




4 


0.012 


0.002 


6 


0.004 


0.002 




6 


0.015 


0.002 


8 


0.005 


0.002 










12 


0.006 


0.003 


4 


2* 


0.057 


0.001 


4* 


0.015 


0.001 




4 


0.048 


0.002 


6 


0.015 


0.002 




6 


0.055 


0.002 


8 


0.017 


0.002 










12 


0.021 


0.002 



The design density is a mixture of two normal distributions. The block size selected most 
often by the Hall, Horowitz and Jing [14] empirical rule is marked by *. 



7. Example. American elk are of high conservation and management 
interest (cf. [1]). This example concerns the relationship between elk and 
landscape characteristics in northern Wisconsin, where elk were reintroduced 
in the early 1990s. There are a total of 514 sampling locations in a 15 km 
X 12 km stand (Figure 4). At each location, a 30 m x 30 m plot was 
surveyed for various plant abundance, including the biomass of grass, forbs 
and shrubs. Among the 514 locations, elk were known to be absent at 250 
plots and present at the remaining 264 plots. One question of interest was to 
compare the biomass indices at locations where elk were absent versus those 
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Table 3 

Estimated nominal coverage of the 90 % GBBB confidence intervals for Po , /3i with sample 
sizes n = 100, 400, 900, scaling factors A„ = 12 with block sizes b„ = 2, 4, 6, and \„ = 24 
with block StZGS On — 4,6,8,12, and a spherical variogram of ranges 2,4,10, based on 
S = 500 simulations, each with M = 1000 bootstrap resamples 











A 7, 


12 








A„ =24 








GBBB 


DSSBB 




GBBB 


DSSBB 


/3o 


/3i 


(3o 


/3i 


/3o 


/3i 


/3o 




100 


2 


2* 


0.876 


0.884 


0.782 


0.82 


4* 


0.88 


0.892 


0.782 


0.84 






4 


0.858 


0.85 


0.788 


0.794 


6 


0.864 


0.888 


0.816 


0.854 






6 


0.774 


0.78 


0.704 


0.702 


8 


0.84 


0.86 


0.794 


0.844 
















12 


0.736 


0.774 


0.686 


0.752 




4 


2* 


0.704 


0.862 


0.58 


0.838 


4* 


0.85 


0.884 


0.72 


0.808 






4 


0.75 


0.832 


0.656 


0.826 


6 


0.838 


0.862 


0.766 


0.826 






6 


0.664 


0.746 


0.556 


0.722 


8 


0.818 


0.846 


0.758 


0.81 
















12 


0.722 


0.756 


0.668 


0.71 


400 


2 


2* 


0.808 


0.896 


0.75 


0.838 


4* 


0.85 


0.902 


0.856 


0.878 






4 


0.814 


0.854 


0.762 


0.812 


6 


0.83 


0.892 


0.85 


0.876 






6 


0.702 


0.746 


0.686 


0.714 


8 


0.818 


0.862 


0.832 


0.832 
















12 


0.724 


0.754 


0.742 


0.728 




4 


2* 


0.648 


0.896 


0.596 


0.86 


4* 


0.822 


0.872 


0.774 


0.866 






4 


0.756 


0.85 


0.726 


0.818 


6 


0.836 


0.846 


0.804 


0.86 






6 


0.67 


0.764 


0.636 


0.716 


8 


0.816 


0.814 


0.804 


0.826 
















12 


0.72 


0.734 


0.68 


0.73 


900 


2 


2* 


0.788 


0.898 


0.802 


0.854 


4+ 


0.856 


0.878 


0.848 


0.892 






4 


0.79 


0.874 


0.81 


0.828 


6 


0.85 


0.858 


0.852 


0.872 






6 


0.688 


0.768 


0.694 


0.75 


8 


0.818 


0.836 


0.816 


0.842 
















12 


0.726 


0.72 


0.724 


0.75 




4 


2* 


0.594 


0.88 


0.566 


0.87 


4* 


0.79 


0.862 


0.758 


0.864 






4 


0.702 


0.82 


0.662 


0.81 


6 


0.806 


0.84 


0.774 


0.84 






6 


0.662 


0.758 


0.584 


0.726 


8 


0.79 


0.824 


0.772 


0.806 
















12 


0.702 


0.74 


0.674 


0.716 



The design density is a uniform distribution. Also reported are the estimated nominal 
coverage of the 90% DSSBB confidence intervals for /3o,/Si- The block size selected most 
often by the Hall, Horowitz and Jing [14] empirical rule is marked by *. 



where elk were present. The answers to the question could help understand 
the forage preferences and the movement patterns of elk. 

We model the biomass indices by a linear regression, z(s) = Pq + Pix{s) + 
e(s), where z{s) is a biomass index, x(s) = l(elk is present at s) and e(-) is 
an error process. Here we account for the spatial dependence in the error 
process and draw inference on (3 by the spatial block bootstrap method 
developed above. The empirical rule shown in Section 6 is used to determine 
the block size. For grass biomass, the block size selected is 6 km. The 90% 
confidence intervals for f3o and /?i are (23.33, 27.04) and (—4.25, —0.35) and 
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Table 4 

Estimated nominal coverage of the 90% GBBB confidence intervals for Po,Pi with sample 
sizes n = 100, 400, 900, scaling factors A„ = 12 with block sizes b„ = 2, 4, 6, and \„ = 24 
with block StZGS On — 4,6,8,12, and a spherical variogram of ranges 2,4,10, based on 
S = 500 simulations, each with M = 1000 bootstrap resamples 









\ 1 r> 

Ati = iJ 






An = 24: 








pi 




po 


pi 


100 


2 


2* 


0.84 


0.888 


4 


0.848 


0.872 






4 


0.816 


0.832 


6* 


0.82 


0.822 






6 


0.74 


0.744 


8 


0.798 


0.79 






— 


— 


— 


12 


0.7 


0.688 




4 


2 


0.648 


0.876 


4 


0.746 


0.858 






4* 


0.676 


0.816 


6* 


0.73 


0.834 






6 


0.562 


0.728 


8 


0.698 


0.81 












12 


0.572 


0.698 


400 


2 


2* 


0.776 


0.876 


4 


0.856 


0.836 






4 


0.76 


0.816 


6* 


0.818 


0.836 






6 


0.66 


0.712 


8 


0.792 


0.812 












12 


0.65 


0.712 




4 


2* 


0.582 


0.898 


4 


0.742 


0.874 






4 


0.66 


0.834 


6* 


0.74 


0.86 






6 


0.538 


0.732 


8 


0.702 


0.822 












12 


0.584 


0.716 


900 


2 


2* 


0.75 


0.878 


4* 


0.83 


0.886 






4 


0.74 


0.828 


6 


0.798 


0.844 






6 


0.616 


0.744 


8 


0.746 


0.796 












12 


0.628 


0.692 




4 


2* 


0.572 


0.872 


4* 


0.744 


0.882 






4 


0.674 


0.804 


6 


0.746 


0.848 






6 


0.592 


0.702 


8 


0.696 


0.822 












12 


0.582 


0.664 



The design density is a mixture of two normal distributions. The block size selected most 
often by the Hall, Horowitz and Jing [14] empirical rule is marked by *. 



the p-value for testing i^o : /3i = is 0.056. There is some evidence that the 
grass biomass was different where elk were present. For the forb biomass, 
the block size selected is 3 km. The 90% confidence intervals for /3o and f3i 
are (14.76, 16.37) and (2.09,4.84) and the p- value for testing i^o : /?i = is 
less than 0.001. There is strong evidence that the forb biomass was different 
where elk were present. For the shrub biomass, the block size selected was 
5 km. The 90% confidence intervals for /3o and (3i are (12.74, 14.77) and 
(0.62,2.98) and the p- value for testing i/o : /3i = is 0.020. There is weak 
evidence that the shrub biomass was different where elk were present. 



26 



S. N. LAHIRI AND J. ZHU 



I 



Ai'ili 



-| 1 1 F 1 r- 

a « I. 4 * A 




4 4 t « 




a r H 4 ■ ^ 



— I — 1 — I — I — I — I — I — 

HI f « ■ i ■ la 





Fig. 2. Boxplots of GBBB variance estimates for combinations of A„ = 12,24, 
n = 100, 400, 900, r = 2,4, and various b„. Two boxplots are placed next to each other 
for /3o and (Si within each subfigure. The sampling design is a uniform distribution. Sim- 
ulation size is S = 500 and resample size is M — 1000. 



If standard linear regression is used [i.e., under the assumption of normal 
i.i.d. e(-)], the p-values for testing i?o : /?! = are 0.030, < 0.0001 and 0.031, 
the 90% confidence intervals for Po are (23.93,26.89), (14.89,16.12) and 
(12.71,14.78), and the 90% confidence intervals for (3i are (-4.81,-0.67), 
(2.21,3.93) and (0.46,3.36), for grass, forb and shrub biomass, respectively. 
After accounting for the spatial dependence, the 90% confidence intervals 
could be either wider or narrower and the p- values could be either above or 
below those from the standard linear regression. 

In the next three sections we give the proofs of Theorem 1 (on the asymp- 
totic distribution of the M-estimators) , Proposition 1 (on inconsistency of 
the DSSBB method) and Theorem 2 (on validity of the GBBB method). 

8. Proof of Theorem 1. Let v(s) = A~^w{s), s € -R„, r]^ = n'U^, C"^ = 
rj^Xn^^'^ , n > 1. Also recall that ttiq^ = sup{||v(s)|| : s € Rn}, n > 1. 
Let tr(i?) denote the trace of a square matrix B. Let denote expectation 

with respect to the joint distribution Px of Xi,X2, Let 

^^=[0,1)*^. Write N = {1,2,...}, ai{t) = f-^i^ and g{t) = t > 0. Set 
A{r,5) = /^°°y(2r-i)d-i^^(y)V2r+<5^y for r G N, 5 G (0, oo), and let m„ = 
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Fig. 3. Boxplots of GBBB variance estimates for combinations of A„ = 12,24, 
n = 100, 400, 900, r = 2,4, and various b„. Two boxplots are placed next to each other 
for j3o and j3\ within each subfigure. The sampling design is a mixture of two normal 
distributions. Simulation size is S = 500 and resample size is M = 1000. 
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Fig. 4. Sampling locations in a study area in northern Wisconsin. The x-axis and the 
y-axis are in the east-west and north-south directions (in m). 
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max{logn,nA~'^}, n>l. For a nonempty set A and a function giA^W, 
let llfl'lloo = sup{|g(x)| : x G yl}. Recall that = cr({Xi, X2, . . .}) and that 
l(-) denotes the indicator function and that C, C(-) denote generic positive 
constants not depending on n and u; € fi. Furthermore, unless otherwise 
specified, limits in order symbols are taken letting n — > 00. 

Lemma 1. Let g:M^M. be a Borel-measurahle function satisfying 
E\g{Z{Q))\ < 00 and Eg{Z{0)) = 0. Also, let ai„(X), i = l,...,n, be 
X -measurable random variables on {Q,T,P) satisfying 



(8.1) ^|ain(X)| = 0(1) asn^oo, a.s. {Px) 

i=l 

and 

n 

(8.2) ^ a^^(X) = 0(1) as 00, a.s. (Px)- 
1=1 

ThenYA=iCLinO^)g{Z{si))^{) in P.\y;^-probability, a.s. {Px.)- 
Proof. This is Lemma 6.2 of [24]. □ 

Proof of Theorem 1. Let j(t) = jQtjj'{Z{si) -'uw(si)'(t - P))du, 
1 <i <n. By Taylor's expansion, the left-hand side of (3.3) equals 



x: w(s,)^(^(s,)) - x: w(s,)w(s,)'(t - /3)ei,(t) 



^.3) 



i=l 



i=l 



= J2w{si)^{Z{s,)) - X] w(s,)w(si)'(t - p)Etl;'{Z{0)) + i?„(t), 
1=1 1=1 

where i?n(t) is defined by subtraction. Note that 
A-ii2„(t) = A-i 



y: w(s,)w(s,)'(t - /3) w'(^(o)) - ^im 

i=l 
n 

K'Y.'^{s^)^{s,)'{k-')'{i;'{Z{s,))-E^|^'{Zm] 



1=1 



+ 



A-iX:w(s,)w(s,)'(A;:i)'{V^'(Z(s,))-4(t)} 



1=1 



AUt-/3) 



i?i„(t) + i?2n(t), say. 
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Next we establish convergence ofn-^J2i=i^i^i)^i^iyWiZ{si))-E'ilj'{Z{0))} 
by using Lemma 1. By Hoeffding's inequality [16] and condition (C.4), it fol- 
lows that for any e > 0, 



Pxl 



n 

n 

1=1 



"'E{ll^(«^)ll'-^x||v(s,)|P} 
1=1 / 

Hence, by the Borel-Cantelli lemma, 

n 

(8.6) n-i5:(||v(s,)f -i^x||v(s,)f) = o(l) a.s. (Px). 

i=l 

By condition (C.l), 

(8.7) £;x||v(si)f = tr(£;xv(si)v(si)') ^tr(i?) asn^oo. 
Hence it follows that 



>e ) <C7exp(-C(e)n/n2'*). 



n 

n 

i=l 



-i^HsOf = 0(1) a.s. (Px) 
and 

n / " \ 

(8.9) n-2^||v(s,)r < n-iE||v(s,)f (n-im2j = o(l) a.s. (Px). 

i=l \ i=l ) 

Thus, by Lemma 1, 



n-i E v(s,)v(s,)'{V'(^(s.)) - E.\^^\Z{m - 

(8.10) 

in P.|x-probability, a.s. (Px)- 

Further, by Proposition 3.1 and Theorem 3.2 of [22] and by a conditional 
version of the Cramer- Wold device, it follows that for n/A^ — > ci G (0, oo]. 



L„^r?J^A-'^/2^v(s,)^(^(s,)) 
i=i 

^.11) 

^ iv(^0, i^a(0)/ci + y" cj(s)Q(s ) ds^ a.s. (Px). 
Also, by arguments similar to those used in (8.5)-(8.7), it follows that 

n 

(8.12) r„ = n-iEv(s0v(s,)'^i7 a.s. (Px). 
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By (8.10)-(8.12), there exist a constant C € (0,oo), a sequence En i and 
sets Bn S n > 1, such that hm„^oo P-\x.{Bn) = 0, a.s. (Px) and on i?^, 

(8.13) \\Lr,\\<Ce~^l\ 

(8.14) ||i?i„(t)|| <nen||A;(t-/3)|| for all t e M*' 
and 

(8.15) \\Tn-H\\<en. 

Without loss of generality, we assume that < log(A„ + 1), n > 1. Hence, 
on the set B^, we may rewrite (8.3) as 

A^t -f3) = {YnX^r^K'^'^Ln + n~\Rin{t) + R2n{t)}]. 

By condition (C.2) and Brouwer's fixed point theorem, this implies that 
P.|x(^n solves (3.3) and |[A;(/3„ - /3)|| < C\-J'^e-^'^) 

(8.16) 

= o(1) a.s. (Px). 
By (8.10) and (8.14) and condition (C.2), it follows that 
(8.17) P.\^mn0n)\\ > A-^/24/2) =o(l) a.s. (Px). 

Hence, by a conditional version of Slutsky's theorem (cf. Lemma 4.1 of [21]), 
it remains to show that 

n 

c-'Y.^{si)^p{Zis,))^N{0,^,) a.s. (Px). 

i=l 

This can be established by using the Cramer-Wold device and by verifying 
the conditions of Theorem 3.2 and Proposition 3.1 of [22]. We omit the 
routine algebraic details to save space. □ 

9. Proof of Proposition 1. For proving Proposition 1, we introduce some 
more notation. Let S'^*(k) be the sum of the DSSBB observations from the 
block i3;*(k), k e /C„, and let sl^\i) denote the sum of the y(sj)'s for 
Sj e I3l^\i), Also, let N** denote the DSSBB sample size, that is, 

N** = \yn{Rn)\, and \eiN = E^N**. Then 

T** = Sl*{\^)/N**. 

kG/C„ 

Since T** is a ratio estimator, we define a linearized version of it, given by 

(9.1) rr= E ^r(k)/iv. 

k6/C„ 
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cr2^Var,(fr) = iV-'A^^,VarJ ^ S^^k)] 

\kG/C„ / 
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(9.2) 



:iV-2A^|/CjVar,(5r(0)) 



Proof of Proposition 1. Part (i) is an immediate consequence of 
Lemma 5.2(ii) of [22]. This fohows by noting that, in this case, the function 
Qi{h.) in the lemma is given by Qi{h) = J /^(s)ds for all h € M"^. To prove 
part (ii), note that for any e > and any 5 € (0, 1), by (9.1), 



<P(|Zj^]|<n5) 



(9.3) 



<P(|Z[f]|<n5) 



+ 



edn 



E 



= /i+/2, say. 
By using the relation E{W) = Ey:_[E.\xW), we have 



E 



(9.4) =E 



n n n 



Y{s,)Yisk)l{sj,Sk G Bl^\i))l{si +Ubn C 

i=ij=ik=i 



n n n 



= E E E ^XCt(s,- - Sfc)l(s„ Sfc G ] (i))l(Si + Ubn C i^n). 
i=l j=lfc=l 

Recall that r„ = Xn/bn- For i j ^ k, by changing variables in the inte- 
grations below, we have 

c„(i,j,k) = ^xfT(sj - Sfe)l(sj-,Sfc G ^|f](i))l(si C Rn) 
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< ^xa(A„(Xi - X2))l(Xi,X2 G X3 + r;,^U) 

= JJ |a(A„(x-y))l(x,yGz + r-i^/)/„(x)/a(y)/a(z)c?xdydz 

(9.5) X faiy + X-h)faiy)fa{z)dsdydz 

<K^ j J J cT(s)l(wGr-iZ^) 

X fa{w + Z + X-'^s)fa(w + z)fa{z)dsdwdz 

X fair'-^u + Z + X~^s)fa{r~^u + z)/a(z) dsdudz. 

By the continuity of /a(-) and the dominated convergence theorem (DCT), 
the right-hand side of (9.5) is asymptotically equivalent to 



(9.6) 



X-dj.-d 



^ j a{s)ds j f^{z)dz = An, say. 
Using similar arguments, it can be shown that 
(9.7) N = n j fl{z)dz{l + o{l)) 

Hence, noting that /C„ = (r„ + l)'^, by (4.6), (9.3)-(9.7), it follows that 



h < 



N-^Xi{rn + 1)" 
e5n 



cn^a{0)+ ^ c„(i,j,k) 



<C 



(9.8) 



Xf , N- ^Xi{r,, + lY ^., 
e5n 



+ 



ieSr' a{s) ds][ I ft{z) dz ) (n^iV-^)(l + o(l)) 



mz)dz mz)dz 



(1 + 0(1)). 



Note that Jf!{z)dz = J ga{^if d^i < 2[(f)3f + b^^ - f)] while 
(//2(z)dz)3 = [/5,(xi)2dxi]3 > [2(|)2i]3 = |^ » //f(z)(iz, as a ^ oo, 
where b = < f for all a > 4. Hence there exist o > 4 and 6 = S{a) € 

(0, 1), r]{a) € (0, 1) (with both 6 and t] sufficiently close to 1) such that with 
e = rja'^, the leading term on the right-hand side of (9.8) is < 1. 
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To obtain an estimate of Ii, define Win = l(AnXj + bnU C Rn), 1 <i <n. 
Then, by the DCT, EWin ^ 1 as n ^ oo. Hence, for any < 5 < 1, by 
Chebyshev's inequahty, 



hm h = hm P(\ll"M < n5) = hm pf V Win < n5 



, i=l 



< hm P 

n— »oo 1 



J2{W^n-EWin) 



i=l 



> n\EWin - 6\ 



< Jirn Var Ij^^in] /[n^EWm - 6f] = 0. 



.1=1 

Thus, it fohows that (4.12) holds with replaced by cj^. To complete the 
proof, it is enough to show that 

This can be shown by using the identity T** = f** + {N - N**)f**/N**, 
the Cauchy-Schwarz inequality and the boundedness of the variables y(-)'s. 
We omit the details. □ 



10. Proof of Theorem 2. 

Proof of Theorem 2. For a positive definite matrix S, let 
denote the probabihty measure corresponding to the A^(0, S) distribution. 
By Taylor's expansion, from (5.11), for any t G M^, we get 

0= [5:(k;t)-c„,(k)] 

kg/C„ 

= [S*n{k-Jn)-Cn{k)] 

ke/c„ 



(10.1) +EE 



w(si)w(sO'(t - /3n)l(si G Bnih; k)) 



x/ ip{Z{si) -uw{siy{t- (3n))du 



J2 [5:(k;/3„)-c„(k)]+A„rX(t-/3„>Xo + A:(t) 

kGK:„ 



say. 



where the remainder term is defined by subtraction. 

Now, by arguments similar to those employed in the proof of Theorem 1, it 
is enough to prove that there exists a nonrandom sequence {en}n>i C (0, 1) 
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with En i as n — > oo such that for any Eq G (0, 1), 



P.|x sup 



(10.2) 
and 

(10.3) 



PJC-I J2 [S:(k;/3„)-c„(k)]Gs) -HB;^c) 



■ o(l), n^oo, a.s. (Px), 



^.|x(n(||C-^P:(t)||>e„||K(t)|| 



for some teW with ||y„(t)|| < | loge„|) > eo) 
= o(l), n— >cx3, a.s. (Px.), 

where we recall that = Xn^^'^r]'^A~^ , and where Vn{t) = Xn^'^A'^{t — (3n)- 
Consider (10.2) first. Let S„, = EkeC„ Var,(C-i5*(k; /?„)) and let (5„ be 
the smallest eigenvalue of S„. By Theorem 17.2 of [3], there exists a constant 
C{p) G (0, oo) such that, on the set > 0}, 



sup 

Sec 



(10.4) 



PACn^ E K(k;/3n)-c„(k)]Gi3) -<I>(S;S. 

<C{p) J2 E4C-^S*^{k-J^)f5-\ 
ke/c„ 



In Lemma 2, we show that 

(10.5) - Sell ^ in P.|x- probability, a.s. (Px)- 
Also, from part (b) of Lemma 2, it follows that 

(10.6) J2 ^*llC~^S:(k;/3„)f ^0 in Pix- probability, a.s. (Px). 

k6/C„ 

Hence, by (10.4)-(10.6), (10.2) follows. 

Next consider (10.3). Note that ^"=1 w(sj)w(si)' = A„r„A^. Hence, 
by condition (C.2), for any t G MP, we have 

IIQ-'K(t)ll 



< 



n 

E E^(«')^(«i)'i(«*^^«(^k;k))-r„ 

kG/C„ i=l 

n 



llK(t)|||xo| 



c 



+ - E E 11^(^0 f+^l(s.Gi?n(/k;k)) 



(10.7) 



kG/C„ i=l 



x(iiK(t)ir+iiK(/3)ir)A-<^^/2||K(t)ii 



+ 
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n-^ E ^v(si)v(sO'{V'(^(si))-Xo}l(si€i?„(/k;k)) 

ke/C„ i=l 
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llK(t)|| 



^RUt)+Rl^{t) + RUt), say. 

Consider -Rt„(t) first. For any e > and A C M*^, let = {x G + 
[—1, Ij'^e C A\. Let ti„ = 6,„/A„, n > 1. Then it is easy to show that uniformly 
in x G -Rq "° , 

(10.8) |{m G Tn : A„x E m + [0, 1)^}| = 6^ + 0(&f i). 

Hence, by (10.8) and the fact that |/Ci„| ~ vol(i?o)^n^n'^ and IT^I ~ vol(i?o)^n' 
we have 

n 

tn = E,n-^ J2 ^ ^n(4;k))l(sj G A„/?o"") 

keK:i„i=i 

n 

(10.9) =n-i5]v(s,)v(s,)'l(s, G A„/?o "")l^in||^nr' E eB„(i;0)) 

= n-i^v(s,)v(s,)'l(s, GA„i?o"") + 
Further, note that by the boundary condition on Rq and (8.6), 



™"'El|v(«jOf E i(s,Gi?n(/k;k)) 

j = l k6^C2n 



(10.10) 

<n~^J2M^j)f\'^^n\\InrX = 0{Un) a.S. (Px)- 

Next, note that by (8.8) and the condition on nion, for any a G (0, oo), 
i5;x||v(si)fl(si^A„i?o"") 

< (Ex||v(si)f+»)l/(2+a)[p^(Xi ^flo"")]'^/(2+a) 

(10.11) 

<m^f+'^)(i?x||v(si)f)-/(2+«).0(n«/(2+-)) 
= o(l) as n — 5- oo. 
Hence by arguments leading to (8.6), a.s. (Px), 



(10.12) 



^Jn-i ^||v(sj)fl(sjGA„Po"") E l(sj GSn(/k;k)) 
I j=i ke/Ci„ 

n 

< ||v(s,-)f l(s,- G A„iio"")|/Ci„||2:„,ri6^ = o(l). 
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Now, by (10.9), (10.10), (10.12) and Chebyshev's inequality, there exists 
ein i such that (10.3) holds with e„ = ei„ and -Rj^(t) replaced by i?J^(t). 
Also, in view of condition (C.2), (10.3) holds for -R2n(t) by Markov's in- 
equality [of the form P=K(|Ty*| > e2n) < -E*|Vr*|/e2n] with e.„ = e2n, where 
£2n = An/^''^^, say. Hence, it remains to prove (10.3) for i?3„(t). To that end, 
let Wji{si) = e'jv{siyv{si)ei['ilj'{Z{si)) - xo], l<j,l<P, i = l,...,n. Then, 
by the independence of the resampled blocks, 



E. 



|x 



Var, E E Wji{si)l{s, G 5„(/k; k)) 



. kG/C„ 



i=l 



(10.13) 



<E.i 



<E.i 



J2 EJn''J2^d^i)^ ^ (4; k) ) 



i=l 



cn~^\inr' E E 



i=l 



2n 



= 0{n-X {<nmny) = o(l) a.s. (Px) 

by Lemma 2(b). Hence (10.3) holds for i?3,„(t). This completes the proof of 
Theorem 2. □ 



Lemma 2. Assume that bn satisfies (5.15). Let g:M^M. be a Borel- 
measurable function satisfying Eg{Z{0)) = and let {z/;„j(X) :1 <i < n}n>i 
be x-fneasurable weight variables with Mn(X.) = max{|t(;„j(X)| : i = 1, . . . , n} < 
oo. Suppose that E\g{Z{0))\'^^'^^ < oo and A{r,6) < oo for some r € N and 
5 G (0, oo). Then 

2r 



(10.14) 



E E ^-ix 

ke/Cn iGXn 



5^z/;„,-(X)5(Z(sj))l(sj e 5„(i;k)) 



= Oi\JCn\\In\bi'[Mn{:^)mnf') a.s. (Px), 
and for any Bn C Rn, 



2r 



^■|x E^™(^)5(^(si))i(sj eS„ 



(10.15) 



.1=1 



= 0(|J„r[M„(X)m„]2n 
«;/iere J„ = {j G Z'^ : P„ n {j + [0, 1)'^} ^ 0}. 



a.s. (Px), 



Proof. Using the exponential inequality in Lemma 5.1 of [22], it can 
be shown that there exists a constant C £ (0, oo) such that 
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n 



Px max Vl(A„X,G{j + [0,l)'^}nS„(i;k)) >Cm„ i.o. 

\lljlli<rf^"/2,ie2:„,kg/c„^ ) 

(10.16) 

= 0. 

Next, grouping the terms in the sum Yl'j=i'^njO^)g{Z{sj))l{sj € i?„(i;k)) 
corresponding to the Sj's in the unit cubes {j + [0,1)*^} PI i?n(i;k)'s and 
applying Theorem 1.4.1.1 of [9] for each i GX„,k€ /Cm one gets (10.14). 
The proof of (10.15) is similar and is omitted. □ 

Lemma 3. Suppose the conditions of Theorem 2 hold. Then 
— Sclj— >0 in P.^x- Probability, a.s. (Px)- 

Proof. For notational simplicity, we prove Lemma 3 assuming that p = 
1. The proof for general p > 1 is similar, except for some additional complex- 
ity in notation. Let 5(1, k) = X;"=i w(sj)^(Z(sj))l(sj G5„(i;k)), ieT„,kG 
/C„, and S„ = Eke^ Jl^n|-i Ei6xJC-i5(i, k))2 - (|T„|-i (C-'5(i, k)))2] . 

Then, by condition (C.2), Lemma 1 and (8.16) and (10.16), it follows that 

(10.17) tn-tn^O in P.|x-probability, a.s. (Px). 

Next, using the regrouping arguments on pages 298-299 of [21], and by 
(10.14) of Lemma 2, it follows that 

^ i?.|x f E [{K'S{i, k)}2 - {E.^x{A-'S{i, k))}']] 

kG/C„ \ ieln ) 

(10.18) 

= 0{\l^n\\ihiml^mih'^) a.s. (Px), 
and by (10.15) of Lemma 2, 

(10.19) 5^ K|x(|X„|-i5^A-i^(i,k)] =0{biml^ml) a.s. (Px). 

kGCn V iGX„ / 

By (10.17)-(10.19), it now remains to show that 

(10.20) lim V |J„ri VP.|x(C^i5(i,k))2 = Se a.s. (Pk). 

Let %n = Ekec,„ l^«r' EiGi„ E.\x{C-^S{\, k))^, j = 1, 2. Then 
-E'x(Si„) = n~^A^|/Ci„,||2'„|~-'- 

Eav,(0)[n^xt'(A„Xi)2 

LieXn 
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+ E - l)^x^i(AnXi)z;(A„X2) 

X c7^(A„(Xi - X2))l(A„Xi, A„X2 € 5„(i;0))} 

= + Si2„,, say, 

where cr^(s) = Eilj{Z{0))ij{Z{s)), s G W^. Note that by (10.8) and condition 
(C.l) [cf. (3.6)], 



^lln = n^r-;— O-^(U) 



(10.21) 



/ ^;(A„x)2/(x) V l(A„xi G i + [0, 1)'^6„) dx 

?;(A„x)2/(x)(ix6^ 



+ 



Ro\Ro ' 



a^(0)i7n~^Af (1 + o(l)) asn^oo. 



1/2 

It can be shown that uniformly in x G and ||s|| < 6n , |{i G 2n : A„x G 

i?„(i;0),AnX + s G i?„(i;0)}| = 6,^(1 + Oibn^^"^)) as n ^ oo. Hence, using 
condition (C.l) [cf. (3.7)], condition (C.5), and change of variables with 
^0 = {x - y : X, y G Rq}, we get 

Si2n = n-U^|/Ci„||J„|-i(n-l) 

(T^(A„u)-t;(A„x)u(A„(x - u)) 



J2 l(A™x, A„(x - u) G Bn{i; 0)) I /(x)/(x - u) dxdu 



In I 



(10.22) 



1n\ J||s||2<fe„ 



/?-""n(flo+A-is) 
x/(x)/(x-A-is) 



v{XnX)v{\nX - s) 



E l(AnX, A„x - s G B„(i; 0)) i dxds 



+ 0(1) 

J a^{s)Q{s) ds + o{l) asn^oo. 



Thus, by (10.21)-(10.22), Ey^tin ^ as n ^ oo. Since |/C2n| = o(|/Ci„|) 
and -Bn(i; k) C Bn{i; 0) for all k G fC2n, a similar set of arguments shows that 
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Ex.'^2n = 0(1) as n — > 00. To complete the proof of (10.20), it is now enough 
to show that 

A„ = tin + S2n - {Ex.tln + ^xS2n) ^ a.S. (Px)- 

This can be done by writing A„ as a U-statistic of order 2 in {Xi, . . . ,X,j} 
and using the arguments developed in the proof of Lemma 5.2 of [22]. We 
omit the routine details. □ 
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